% Table 9 in "Commodity Prices, Convenience Yields, and Inflation"
% March 2011

clear all;

load cyj.dat;                   % convenience yields on 23 commodities
load us_intrate_monthly.dat;    % US 3-month T-bill rate
load us_cpi_monthly.dat;        % US CPI indeces (all items, lfe etc.)
load imf_index_monthly.dat;     % IMF aggregate commodity price index
load usd_monthly_ave.dat;       % USD trade-weighted index

i=us_intrate_monthly;
P=us_cpi_monthly;
SA=imf_index_monthly;
er=usd_monthly_ave;
nn=rows(P);

rhat=2;         % number of principal components
nwl=4;          % number of lags in Newey-West HAC estimator
maxp=4;         % maximum lag for AIC
hh=[1 3 6 12];  % forecast horizon

for hi=1:4;
h=hh(hi);

TT=177;
act=[]; forc1=[]; forc2=[]; forc3=[];
mse1=[]; mse2=[]; mse3=[];
for j=1:nn-TT-h;
    ds_h=100*(log(SA(1+maxp+h:TT+j,1))-log(SA(1+maxp:TT+j-h,1)));
    [ehat,Fhat,lamhat,ve]=pc(standard(cyj(1:TT+j,:)),rhat);
    const=ones(rows(ds_h),1);
    nr=rows(const);
    aic1=zeros(maxp,1); aic2=zeros(maxp,1); aic3=zeros(maxp,1);
    ds_e=zeros(nr,maxp); p1_cy=zeros(nr,maxp); p2_cy=zeros(nr,maxp);
    dx=zeros(nr,maxp); ir=zeros(nr,maxp); ds_f=zeros(1,maxp); 
    p1_cyf=zeros(1,maxp); p2_cyf=zeros(1,maxp);
    dx_f=zeros(1,maxp); ir_f=zeros(1,maxp);
    for ip=1:maxp;
        ds_e(:,ip)=100*(log(SA(2+maxp-ip:TT+j-h+1-ip))-log(SA(1+maxp-ip:TT+j-h-ip)));
        p1_cy(:,ip)=Fhat(2+maxp-ip:TT+j-h+1-ip,1);
        p2_cy(:,ip)=Fhat(2+maxp-ip:TT+j-h+1-ip,2);
        dx(:,ip)=100*(log(er(2+maxp-ip:TT+j-h+1-ip))-log(er(1+maxp-ip:TT+j-h-ip)));
        ir(:,ip)=i(2+maxp-ip:TT+j-h+1-ip,1);
        ds_f(:,ip)=100*(log(SA(TT+j+1-ip))-log(SA(TT+j-ip)));
        dx_f(:,ip)=100*(log(er(TT+j+1-ip))-log(er(TT+j-ip)));
        ir_f(:,ip)=i(TT+j+1-ip,1);
        p1_cyf(:,ip)=Fhat(TT+j+1-ip,1);
        p2_cyf(:,ip)=Fhat(TT+j+1-ip,2);
        res=nwest(ds_h,[const p1_cy(:,1:ip) p2_cy(:,1:ip) ds_e(:,1:ip)],nwl);
        res2=nwest(ds_h,[const dx(:,1:ip) ir(:,1:ip) ds_e(:,1:ip)],nwl);
        res3=nwest(ds_h,[const ds_e(:,1:ip)],nwl);
        aic1(ip)=log(res.sige)+2*(rows(res.beta)-1)./rows(const);
        aic2(ip)=log(res2.sige)+2*(rows(res2.beta)-1)./rows(const);
        aic3(ip)=log(res3.sige)+2*rows((res3.beta)-1)./rows(const);
        indep1=[1 p1_cyf(:,1:ip) p2_cyf(:,1:ip) ds_f(:,1:ip)];
        indep2=[1 dx_f(:,1:ip) ir_f(:,1:ip) ds_f(:,1:ip)];
        indep3=[1 ds_f(:,1:ip)];
        forc1(j,ip)=indep1*res.beta;
        forc2(j,ip)=indep2*res2.beta;
        forc3(j,ip)=indep3*res3.beta;
    end;
    [aicmin1,minind1] = min(aic1);
    [aicmin2,minind2] = min(aic2);
    [aicmin3,minind3] = min(aic3);
    forc4(j,:)=100*(log(SA(TT+j))-log(SA(TT+j-h)));
    act(j,:)=100*(log(SA(TT+j+h,1))-log(SA(TT+j,1)));
    mse1(j,:)=(forc1(j,minind1)-act(j))^2;
    mse2(j,:)=(forc2(j,minind2)-act(j))^2;
    mse3(j,:)=(forc3(j,minind3)-act(j))^2;
    mse4(j,:)=(forc4(j)-act(j))^2; 
    mae1(j,:)=abs(forc1(j,minind1)-act(j));
    mae2(j,:)=abs(forc2(j,minind2)-act(j));
    mae3(j,:)=abs(forc3(j,minind3)-act(j));
    mae4(j,:)=abs(forc4(j)-act(j));
end;

bench1=sqrt(mean(mse3));
bench2=mean(mae3);
fprintf('    h = %6.4f\n',h);
fprintf('    RMSE\n')
disp([sqrt(mean(mse1))./bench1 sqrt(mean(mse2))./bench1 sqrt(mean(mse4))./bench1])
fprintf('    MAE\n')
disp([mean(mae1)./bench2 mean(mae2)./bench2 mean(mae4)./bench2])

end

